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Abstract. - A Brillouin-Wigner perturbation theory is developed for open electromagnetic sys- 
tems which are characterised by discrete resonant states with complex eigenenergies. Since these 
states are exponentially growing at large distances, a modified normalisation is introduced that 
allows a simple spectral representation of the Green's function. The perturbed modes are found by 
solving a linear eigenvalue problem in matrix form. The method is illustrated on exactly solvable 
one- and three-dimensional examples being, respectively, a dielectric slab and a microsphere. 
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^ r*| . Eigenvalue problems are often not analytically solvable, 
Q-and approximate solution schemes are employed, which 
use a known (ideally analytic) solution of a similar prob- 
„ lem, differing from the problem of interest by the 'per- 
^turbation'. The most well known scheme is the Rayleigh- 
t" Schrodinger perturbation theory pQ , which writes the so- 
£^ lution as series expansion in terms of the perturbation in- 
volving sums over all eigenstates, and thus introduces the 
•notion of perturbation orders. However, the convergence 
properties of such series limit the treatable perturbation 
strength, and the mathematical complexity limits the ex- 
^Hpansion typically to a few orders. One way to overcome 
^this limitation is to truncate the infinite complete set of 
• i— i unperturbed eigenstates to a finite number and then use 
^matrix diagonalisation algorithms to treat the perturba- 
^tion within this basis exactly, i.e. to all orders. In quan- 
tum mechanics, this method is called Brillouin-Wigner 
(BW) perturbation theory [2]. If the unperturbed spec- 
trum contains bound states only (as for the harmonic os- 
cillator), it allows to treat the effect of perturbation to any 
required accuracy by taking a sufficiently large number of 
states. However, if a continuum of scattering states is im- 
portant, the use of the BW theory becomes problematic. 
In electromagnetic systems, such a continuum is often the 
dominating part of the spectrum and discrete bound states 
(evanescent waves) are the exception. 

Both Ray leigh- Schrodinger and BW perturbation theo- 
ries are well established when treating Hermitian prob- 



lems which describe conservative systems having real- 
valued eigenenergies and orthogonal and normalised wave 
functions. In open systems instead, an eigenstate ac- 
quires a certain linewidth which characterises its finite life- 
time, so that the eigenenergy becomes complex- valued [3]. 
The corresponding eigenfunctions are not orthogonal and 
grow exponentially in the outer space so that they cannot 
be normalised by integrating their square modulus @HB]. 
Both the energies and the wave vectors of these states 
have negative imaginary parts resulting in a wave func- 
tion oc exp[— T(t — r/v)] (T = — Imcj, v is the phase 
velocity). The latter expression provides a physical pic- 
ture of the exponential growth as a wave front propagat- 
ing away from the system which has been excited at an 
earlier time. In this way stationary bound states become 
quasi-stationary [5] in open systems due to leakage to the 
outside. In the literature these states are also known as de- 
caying states [7], resonant states (RS) [8], leaky modes [9], 
quasi-guided modes [10], quasi- normal modes [TTj [T2], etc. 
Moreover, going away from the real-energy axis into the 
complex energy plane, a continuum of stationary scatter- 
ing states, which is problematic in perturbation theory, 
can be effectively eliminated from the spectrum and re- 
placed by a countable number of RS [6] . Formally, this is 
achieved by imposing boundary conditions of no incoming 
waves, i.e. no waves travelling towards the system. 

The Ray leigh- Schrodinger perturbation theory of RS 
has been developed for quantum-mechanical [7] and elec- 
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tromagnetic systems [12]. For weak perturbations, an ap- 
proximate BW theory limited to degenerate modes only 
has been reported [13]. Also, an exact method of cal- 
culation of RS perturbed by a point scatter has recently 
been developed [TJ. However, sharp resonances in three- 
dimensional structures cannot be calculated by known 
methods in electrodynamics. Thus it is especially impor- 
tant to develop a BW perturbation theory, which could 
treat such resonant states under a perturbation of arbi- 
trary strength and shape. Such a theory was only partially 
formulated in electrodynamics more than 40 years ago [6] 
and later on proposed also in the quantum-mechanical 
context p~5l[l6], but has not been completed neither ap- 
plied to numerically solve relevant problems. 

In this Letter we present a BW perturbation theory of 
open electromagnetic systems and validate it on exactly 
solvable one-dimensional (ID) and three-dimensional (3D) 
examples, for both weak and strong perturbations. Al- 
though the RS have complex energies and wave functions 
exponentially increasing outside the system [5], orthogo- 
nality and normalisation can be established. Then, the 
expansion of the perturbed system into the RS of unper- 
turbed system can be reduced to the standard BW form 
of diagonalisation of the sum of a diagonal matrix of the 
unperturbed system and a perturbation matrix. 

To make the idea clear and the derivation transparent 
we present our BW theory first on an effective ID prob- 
lem taking a homogeneous dielectric slab as unperturbed 
system and then, after bringing it to a general form, use 
it for a dielectric microsphere. 

We consider a dielectric slab in vacuum with thickness 
2a which is described by a real dielectric constant along 
the spatial coordinate z [see the inset in fig. 1(a)]: 



e(z) 



1, 



\z\ < a, 
\z\>a. 



(1) 



The permeability is /x = 1 everywhere throughout this 
work. For transversal eigenmodes E n = £E n /y/S (S is 
the normalisation area in the x?/-plane) having zero in- 
plane wave number, k x = k y = 0, and time-dependent 
part given by ex.-p(—ick n t) with the complex frequency 
cfc n , the Maxwell equation is reduced to 



E n {z) = 



(2) 



with the electric field E n changing smoothly {i.e. having 
a continuous derivative) across the boundaries at z = d=a. 
Assuming the boundary conditions (BC) of only outgoing 
waves results in the eigenfunctions 

{-l) n A n e- ik " z , z < -a 

E n (z) = { B n [e^knz + (_^n e -v^n*] 5 \ z \<a 

z > a 

(3) 



and the wave number eigenvalues 



1 



2a^/T s 



(yrn-Hna), n = 0, ±1 , ±2 . . . , (4) 



where a = (>Jel + l)/{^/el — 1). The eigenfunctions E n 
are orthogonal and normalised in a well-known way [6] 
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E n (a)Em(a) + E n (-a)E m (-a) 



i{k n ~\~ km) 
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-ik n a 
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Let us also construct the unperturbed Green's function 
(GF) G(z, z'\ k) which satisfies the equation 
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G(z,z';k) = 5(z-z'), 



(7) 



and obeys the same BC of no incoming waves. As in 
the case of a square-barrier potential in quantum mechan- 
ics [T7] , this GF can be calculated analytically. Inside the 
slab (— a < z, z' < a) it has the following explicit form: 



G(z, z'\ k) 



-|_ e ig(a-2>) _|_ ae -iq{a-z > ) 

2iq e^ a + ae~^ a 

x- . . (8) 



where z > = max(z,z / ), z < = min(^,z / ), and q = k-^/e^. 
The analytic continuation of the GF into the complex k- 
plane has simple poles at k = k n and k = 0. Then, similar 
to the resolvent expansion of a Hermitian operator in the 
energy plane [19], it can be written, using the Mittag- 
Leffler theorem, as a sum of partial fractions [7([T8]: 



G(z,z';k) = 



E n (z)E n (z') 



2kyi(k kfi) 2zk 



- = T 



E n (z)E n (z') 
2k(^k kfi) 



(9) 

where the residua were expressed in terms of the 
normalised eigenmodes eq.® and a sum rule 
J2 n E n (z)E n (z / )/k n = i following from the high- 
frequency asymptotics of the GF [20] was used. 

Let us now consider an arbitrary real perturbation 
Ae(z) of the dielectric constant inside the slab. The wave 
equation for the perturbed resonant states £ v and wave 
numbers k v is given as in eq. (|2]) by 



— +[s(z) + As(z)]>Z\£„(z) = 



(10) 



with BC £ v (z) oc exp(zx^|z|) at \z\ > a. The new GF 
G(z, z'\ k) is related to the unperturbed one via the Dyson 
equation 

g(z,z , \k) = G(z,z , ]k) 

-k 2 ^ G(z,z"]k)Ae(z")g(z",z']k)dz". (11) 

J — a 

On the other hand, a spectral representation of the per- 
turbed GF, similar to eq. (j9]), can be used. It is derived 
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in the Appendix for an arbitrary profile of the dielectric 
function inside the slab. Then, following ref. [7] we equate 
residua at the poles k = k v in eq. (fTT|) . This results in the 
following relationship between the unperturbed and per- 
turbed modes: 



E n (z)C a E n (z')Ae(z>)£„(z')dz> 
2{k n /x v - 1) 



(12) 



In the interior region \z\ < a which contains the perturba- 
tion, the perturbed RS can be expanded into the unper- 
turbed ones, exploiting the completeness of the latter, 



e v {z) 



b ni/ E n (z) . 



(13) 
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Substituting this expansion into eq. ([12]) and equating co- 
efficients at the same basis functions E n (z) results in the 
following matrix equation: 



2b n 



k n / x v 



nm u mv 



(14) 



Introducing new coefficients c nv = b nv \ 
the matrix problem, we arrive at 
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r h n to symmetrise 



(15) 



Equation ([T5]) is the central result in our BW perturba- 
tion theory of resonant states. Similar to the BW the- 
ory of bound states in quantum mechanics, eq. (fT5|) rep- 
resents a linear eigenvalue problem which requires diago- 
nalisation of a matrix depending only on the unperturbed 
spectrum and the perturbation [21]. However, instead of 
real eigenenergies and a Hermitian perturbation matrix 
characteristic for conservative systems, we are dealing here 
with complex- valued bare and perturbed eigenvalues, l/k n 
and 1 / k v , respectively, and a complex symmetric pertur- 
bation matrix V nm / (2\/k n k m ), where 



V nm = J A£(r)E n (r)-E m (r) dr. 



(16) 



While from the formulation of the method both the unper- 
turbed system and the perturbation Ae(r) can be located 
inside any finite area, we found that good convergency is 
achieved for a perturbation within the area of the unper- 
turbed system. 

To illustrate the accuracy of the method, we use as un- 
perturbed system the slab with ^/ej = 1.5. The unper- 
turbed discrete spectrum is given by eq. (|4]) and shown in 
fig. l(a,b) by circles with a dot. We introduce a strong 
non-symmetric perturbation As =10 and 100 in the layer 
(a/2 < z < a) [see the inset in fig. 1(a)], thus increas- 
ing in the latter case the refractive index by an order of 
magnitude. The exact eigenvalues x® xact are given in the 
Appendix by eq. (A. 9). They are compared in fig. 1 (a,b) 
(squares) with our numerical result (crosses) calculated 



Fig. 1: (colour online) Resonant states of a dielectric slab. 
Exact and calculated via eq. (p"5)) wave numbers x u (a,b) along 
with the relative error of the calculation (c,d), for different 
truncation N of the matrix problem. Unperturbed: RS of a 
dielectric slab in vacuum (black circles with a dot). The slab 
with s s = 2.25 occupies the area —a < z < a. Perturbed: 
Exact (red rectangles) and calculated (blue crosses) RS wave 
numbers x v of the same slab but with a higher refractive index 
in the layer (a/2 < z < a), the change there being As — 
10 (a,c) or As = 100 (b,d). Inset: profiles of the dielectric 
constant of the unperturbed and perturbed systems. 



according to eq. (fT5j) where the sums have been truncated 
at —N < n < N, yielding a matrix dimension of (2Af + 1). 
The relative error \k v ) x® xact — 1| is quantified in fig. 1 (c,d) 
for different N. States close to the truncation border 
(i.e. with large values of are always poorly calculated. 
However, as the matrix size increases, a state with a given 
v is moving closer to the middle of the matrix, so that its 
relative error decreases as N~ 3 . The exact eigenfunctions 
(the electric field £ v ) are reproduced to a similar accuracy. 

Writing the perturbation in the form of eq. (fTB]h the 
bare eigenmodes E n (r) are supposed to be normalised as 
given by eq. (|5j) in the ID case. In quantum- mechanical 
systems, normalisation of RS was discussed e.g. in 
refs. [7J[8l[22]. A similar approach was used to normalise 
eigenmodes in spherically symmetric open electromagnetic 
systems [6|ll2pi3j . Some important details of such normal- 
isation are given in the Appendix along with the spectral 
representation of the Green's function, see eqs. (A. 23) and 
(A. 24). Below we show, however, how this normalisation 
can be done in the general 3D case. 

The electric field modes E n obey Maxwell's equation 



Vx VxE n (r) = fc£e(r)E n (r), 



(17) 



where the dielectric constant e(r) is assumed to deviate 
from unity only inside a sphere of radius a. The perturba- 
tion Ae(r) is located within the same sphere. Multiplying 
eq. (|T7|) by E m (r) and integrating the dot product over the 
sphere of radius R ^ a by parts, we obtain the orthogo- 
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Fig. 2: (colour online) Resonant states of a sphere, with the 
angular momentum I = 20 (TE modes). Unperturbed: Exact 
RS wave numbers k n of a homogeneous dielectric microsphere 
in vacuum, with radius a and refractive index 1.5 (black circles 
with a dot). Perturbed: Wave numbers x v of a microsphere 
with the same radius a but increased dielectric constant (Ae = 
5: red rectangles and blue crosses +; As = 25: red circles 
and blue crosses x). Results of eg. ([T5]) are compared with 
exact values. Left inset: profiles of the dielectric constant of 
the unperturbed and perturbed systems. Right inset: relative 
errors of calculation with As = 5 for different truncation N of 
the matrix problem. 



nality condition 

= (k 2 n - k 2 m ) f dve{v)E n {v) • E m (r) 
Jr 

dE n 



dS[E n 



dE m 

dr 



dr 



(18) 



where the last integral is taken over the surface of the R- 
sphere. Both integrals in eq. (fT8|) are independent and 
divergent at R — >• oo. There is, however, an exact cancel- 
lation of these ^-dependencies thus removing the diver- 
gencies [23]. In a degenerate case k n = k m (n / m), both 
the volume and the surface integrals vanish by symmetry 
showing the orthogonality. 

Equation (fT8|) is now used as a starting point for the 
eigenstate normalisation. In the outside range r > a, 
where e(r) = 1, we take advantage of the functional be- 
haviour E n (r) = F n (k n r) with an auxiliary function F n 
being a superposition of vector spherical harmonics [24] . 
The analytic continuation E(fc,r) of the eigenfunction 
E n (r) around the point k n in the complex /c-plane is intro- 
duced in such a way that E(fc,r) = F n (fcr) and satisfies 
eq. (fT7]) (with k n being replaced by the k). The normali- 



sation of E n (r) is then defined as follows: 



dve(v)E 2 n (v) 



lim 

k->k n 



V dr dr_ 

rh ri rh 



■ (19) 



It is again i?-independent, so that eq. (fT^jl can be evalu- 
ated for R = a. Furthermore, due to the above mentioned 
functional behaviour, the limit k — > k n can be taken ex- 
plicitly reducing the last term in eq. (|19p to 



1 

2fcf 



dS 



dr 



rE n 



d 2 E n 

dr 2 



dE n 

dr 



(20) 

To illustrate the 3D case, we take as unperturbed sys- 
tem a dielectric microsphere of radius a and refractive in- 
dex equal to 1.5. We consider a perturbation (sketched 
in the left inset of fig. 2) which uniformly increases the 
refractive index inside the sphere (As = 5 and 25) leav- 
ing its radius unchanged. Since both the unperturbed and 
perturbed systems have spherical symmetry, the angular 
momentum I is a good quantum number and can be used 
to classify the RS. Additionally, eigenmodes can be clas- 
sified as transverse electric (TE) and transverse magnetic 
(TM). These TE and TM modes are not mixed by our 
perturbation. The dielectric microsphere has a rich spec- 
trum containing both leaky and high-Q whispering gallery 
(WG) modes, the latter having a small imaginary part of 
h v . The difference between these two types of modes be- 
comes more prominent at larger values of /, see the spec- 
trum for I = 20 in fig. 2. For TE modes, the exact solutions 
(squares and circles) are well reproduced by our method 
(crosses). To calculate the eigenstates shown in fig. 2, the 
matrix problem was truncated to TV = 200. The accuracy 
quickly improves as the matrix dimension (27V) increases, 
again with a relative error oc TV -3 , as is seen in the right 
inset of fig. 2. 

In conclusion, a Brillouin-Wigner perturbation theory of 
resonant states in open electromagnetic systems has been 
developed and validated on exactly solvable ID and 3D 
examples. Although the wave functions are exponentially 
growing at large distances and thus cannot be normalised 
in the usual way, the eigenvalue problem is reduced to 
the standard linear form, typical for conservative systems. 
Practically, it only requires diagonalisation of a complex 
symmetric matrix which consists of the bare spectrum and 
the perturbation. The simple spectral representation of 
the Green's function in terms of the calculated perturbed 
resonant states allows an easy access to measurable quan- 
tities such as emission, scattering, and transmission. 



The authors thank P. Borri, A. L. Ivanov and L. 
Chantada for valuable discussions. E.A.M. acknowledges 



p-4 



Brillouin-Wigner perturbation theory in open electromagnetic systems 



financial support of WIMCS. W.L. acknowledges support 
from BBSRC, grant BB/E005624. 

Appendix 

Spectral representation of the Green's function 

in ID. — The central quantity for calculating e.g. trans- 
mission and emission properties of a dielectric system is 
the corresponding GF which in ID satisfies eq. ([7]) with 
an arbitrary profile e(z) of a dielectric slab embedded in 
vacuum [i.e. e(z) = 1 for \z\ > a]. It can be solved in the 
form 

E L ( Z< ,k)E R ( Z> ,k) 



G{z,z'\k) = 



(A.l) 



W L , R (k) 

with EL,(z,k) and E R (z,k) being the left- and right- 
exponential homogeneous solutions, 

E L (z,k) = e~ ikz for z < -a , 

E R (z, k) = e +ikz for z > -fa , (A.2) 

so that the GF obeys the BC of outgoing waves at either 
side, for real positive wave number k. The Wronskian 
W LiR (k) = E L (z, k)E' R (z, k) - E' L (z, k)E R (z, k) in the de- 
nominator takes care of the delta inhomogeneity in eq. ([7]) • 
The analytic continuation of the GF into the com- 
plex /c-plane has simple poles at k n where the Wron- 
skian vanishes, i. e. where left- and right-exponential so- 
lutions agree (up to a factor). This defines the eigenmodes 
E n (z) oc £x(z,fc n ) oc E R (z,k n ). They have an exponen- 
tial behavior on both sides of the slab, 



E n (z) oc ex.p(ik n 



for I z I > a . 



(A.3) 



For any real dielectric function, the eigenstates come in 
pairs [E n {z) and E*(z)], giving rise to wave numbers 
which are positioned symmetrically with respect to the 
imaginary k axis. The states with Im k n < which we 
call RS increase exponentially outside the slab and there- 
fore the orthogonality relation and normalisation have to 
be modified. 

The orthogonality is established by integrating 
Maxwell's equation (|2]) for E n (z), multiplied with E m (z), 
over a finite interval (— R, +R) with R > a. Integrating 
by parts the derivative term, and subtracting the inter- 
changed expression n o m gives 

i(k n - k m ) [E n E m ] s + (k 2 n - k m )N nm = (A.4) 

with the overlap integral 

r+R 

N nm = / dze(z)E n (z)E rn (z) (A.5) 

J-R 

and a 'surface' term, [F] s = F(R) + F(-R). Note that 
the independence of the surface term [E n E m ] s and that 
of the overlap integral N nm cancel each other exactly. We 
choose normalisation of the RS in such a way that this 
property is preserved also for n = m: 



[E n E. 



™>\S 



i{k n + km) 



(A.6) 



so that the RS do not depend on R. The normalisa- 
tion constant is chosen such that for bound states (having 
Rek n = and Imfc n > 0) the standard normalisation 
N nn = 1 is restored, as here the surface term can be made 
small for large R. A careful investigation shows that the 
Wronskian derivative with respect to k equals 2fc n , so that 
the GF takes the form of eq. ©. 

Eigenmodes of the perturbed slab in fig. 1. — 

The profile of the dielectric constant of the perturbed slab 
sketched in the inset to fig. 1 has the following form 



e(z) + Ae(z) 



1. 



1, 



z < — a , 
-a < z < b , 
b < z < a , 
z > a . 



(A.7) 



The eigenmodes are found as solutions of the Maxwell 
equation ([TO]) with the BC of outgoing waves and have 
the following explicit form 



£u(z) 



Fe ix v z 



z < —a , 

. Ce -iy/e7* v z ^ - a <z<b, 
- Ee-ty/^""* , b<z<a, 
z > a . 
(A.8) 

The six coefficients in eq. (A.8) are found (for each v) us- 
ing the continuity of the E v (z) and £' v (z) on the three 
boundaries (z = —a, z = 6, and z = a) and the nor- 
malisation condition eq. (A.6). The set of six equations 
following from the above mentioned BC not only estab- 
lishes the relations between the coefficients but also leads 
to the following secular equation for the wave numbers k v : 



aPf{xu)g{xu) - 1 = 



a/3 - 1 



Pg(x v ) - af(x v ) , (A.9) 



where 



a 



(3 



(A.10) 



and the functions f(k) and g(k) are defined as follows: 

f(k) = exp [-2iy/T 8 k{a + b)} , (A.ll) 

g(k) = exp [-2iy/e^k(a - b)] . (A.12) 

Derivation for the 3D case with spherical sym- 
metry. — We consider a dielectric sphere of radius a with 
a dielectric constant of spherical symmetry, e(r), embed- 
ded into vacuum. The solutions of Maxwell's equation 
can be classified into TE and TM modes, and the angular 
momentum is a conserved quantity. The Green's function 
can be separated into /-diagonal contributions, 

1 00 2/ -1-1 

G(r,T>;k) = — J2Gi(ry ; k)——P l (cose) (A.13) 



1=0 



where O denotes the angle between r and r 7 , and Pi(x) 
are the Legendre polynomials. For simplicity, we consider 
only the TE case below. The GF obeys the equation 

d? _ Z(Z + 1) + k2£{r) \ Gi ^ r ,. k) = 5 ( r _ r ,^ (A M) 
ar z r z 
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which differs from the ID version eq. ([7]) only by the cen- 
trifugal term oc 1/r 2 . For the homogeneous equation, we 
define a solution /fc(r) which is regular at the origin, 

fk(r) = r l+1 for r ^ (A.15) 

and another one, #fc(r), which has outgoing wave BC: 



9k(r) 



-\-ikr 



for r 



oo . 



(A.16) 



As in the ID case, the GF can be constructed using both 
solutions, 

G l (ry-,k)= fk £<l 9kir > ) . (A.17) 
W{fk,gk} 

The analytic continuation into the complex k plane has 
poles wherever the Wronski determinant has zeroes which 
defines the wave numbers of the resonant modes, 



W{f k ,gk} = at k = k n . 



(A.18) 



At k n , both solutions fk and g k agree up to a factor. In- 
tegrating Gi(r,r';k')/(k' — k) over a large circle in the 
complex k' plane gives zero, but can be decomposed into 
pole contributions with the result 



G l (r,r';k) = J2 



fkn ( r ) 9k n (r') 
(k - k n ) w n 



(A.19) 



(note that in 3D, there is no trivial pole at k = 0). The 
Wronskian derivative, w n = dW '/ 'dk\k= kn , is calculated 
integrating Maxwell's equation by parts over the finite in- 
terval (0, R) with R> a which gives 

dW{f k: g k } f (p , dg' k (R) fl(m dg k {R) 
~ Jk{R) — Jk{ H ) 



dk 



dk 



dk 



+2k dre(r)g k (r)f k (r), (A.20) 
Jo 

owing to the /c-independent regular solution eq. (A.15) at 
r — )• 0. A further simplification is possible due to R > a 
where e(r) = 1 holds. Here, the outgoing wave solution 
depends on r and k only in product form, g k (r) = g(x = 
kr). Then, the first term on the r.h.s. of eq. (A.20) can be 
rewritten as 



f k (R) [g'(x) + xg"(x)} - f' k (R)Rg'(x)\ x __ 



kR 



(A.21) 



At the poles, we have fk n ( r ) — c n • gk n (r), and can take 
9n(r) =gk n (r) everywhere, 



w n = g{x)g\x) + x [g{x)g"{x) - {g\x)f 



-k n R 



+2k n I dre(r)gl(r) 
Jo 



in conjunction with 



G,(r,r^) = ^ 



9n{r) g n {r') 
(k - k n )w n 



(A.22) 



(A.23) 



Equation (A.22) can be considered as the proper normal- 
isation of the resonant modes. Its l.h.s. is formally equiv- 
alent to their orthogonality (valid for k n ^ km), 







g n (R)g' m (R) ~ 9n(R)9m(R) 
km k n 



+ (k n + k m ) I drs(r)g n (r)g m (r), (A.24) 
Jo 

if the limit k m — >• k n is performed. 
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